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1. Introduction 

The spectrum of cosmic rays (CR) in the local interstellar medium (within about 1 kpc from 
the Sun) is of interest as a complement to direct measurements and as a probe of solar modulation 
which affects particles below a few GeV. y-rays from CR protons and heavier nuclei interacting 
with interstellar gas are an ideal probe of local CR; bremsstrahlung from CR electrons and positrons 
is also important at low energies and must be accounted for too. In (1; 2) we presented preliminary 
results based on earlier emissivity data (3). (4) addresses other aspects of this topic, also using 
those emissivities. 

Recently a new and precise determination of the local y-ray emissivity using Fermi -LAT data 
has been made (5), studying regions at Galactic latitudes out of the plane; in particular the contri¬ 
butions from atomic, molecular and ionized hydrogen were separated, and the emissivity of atomic 
gas traced by the 21-cm line can be used as the most reliable data for analysis. (5) also provided 
an extensive analysis of the emissivities to derive the interstellar proton spectum. Here we pursue 
this project with the new emissivities, with various innovations including new cross-sections and 
analysis techniques. 

2. Emissivity matrix representation 

2.1 Response Matrices 

The observed emissivity is the sum of hadronic (pion-production) and leptonic (bremsstrahlung) 
contributions: q tot (Ey) = qi ia d{Fy ) + qi ep (Ey). The gamma-ray emissivity predicted for the input 
CR spectra can be written in matrix form E^ red = 'L process LjM^1)(d). where £k is the emissivity in 
the k’th energy bin, and Ij(9) is the j'th momentum sample point of the CR spectrum model with 
parameters 6. The processes run over hadronic and leptonic components, with I, (6) including CR 
protons, helium, electrons and positrons. Mj/ C can include the dispersion in the Fenni-LAT y-ray 
energy measurements. The matrix Mjk is pre-computed for a given binning of data and model, so 
that the response can be computed very fast for a large number of model parameters. The method is 
general and can be extended to include cosmic-ray direct measurements, synchrotron emissivities 
etc. Once the matrices and data are defined, the problem is a purely mathematical one. 

2.2 Energy dispersion 

The Fenni-LAT energy measurement has a precision of about 15%, becoming worse at low 
energies. Accounting for the energy dispersion is especially important below 200 MeV, where pho¬ 
tons are both lost from the measured range and gained from higher and lower energies, depending 
on the input spectrum. The energy dispersion is defined as p(E meas \E true \ where E true ,E meas are 
the true and measured y-ray energies respectively. p(E meas \E true ) is obtained in matrix form using 
the Fermi -LAT Science Tools, for the same event selection and response as used to derive the emis¬ 
sivities. The approximation is made that the response averaged over the full sky is appropriate; the 
variations are small and average out making this a good approximation. The measured spectrum 
is also weighted by the Fermi -LAT effective area which is a function of energy, varying rapidly 
at low energies. The full transfer function from the model proton and helium spectra to measured 
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emissivities in energy bands can be expressed as a single matrix, allowing fast computation for 
statistical analysis. 

3. Bayesian analysis 

The likelihood function is L{data\Q) = exp(—£,t (e£ red ^ — ef >s ) 2 /2o[) where £% bs is the ob¬ 
served emissivity and Cj. is the error estimate for the k’th energy bin. By Bayes theorem the 
posterior probability of the model is P(9\data) = - where P(9) is the prior probabil¬ 
ity for the parameters 6 defining the model, and P(data) is a normalizing factor (known as the 
evidence ) such that fP(0\data)d N 6 = 1. Given the joint probability distribution of the param¬ 
eters, we can marginalize to obtain the distribution of any subset, including single parameters: 
P(9i) = Jo/e, P(9\data)d N 1 0. Mean and standard deviations are 0, = Jo/e, 0jP(6\data)d N 1 0, 
o(9j) = \/( Jo/e,^' — 0j) 2 P(0\data)d N ~ l 6). Since we are actually more interested in the CR 
spectrum than in a particular set of parameters, which are anyway highly correlated, it is useful 
also to compute the mean and standard deviation of the gridded CR spectra tij = n(pj): 

hj = Jnj(9)P(9\data)d N 9, a(nj) = \/(f(n/(d) — hj) 2 P{9\data)d N 9). This defines a band 
containing the range of CR spectra. This method is therefore a deconvolution using a basis of 
parameterized spectra 1 . The range of emissivities corresponding to this range can then be computed 
(for each process) for comparison with the observed values, to illustrate the uncertainties. 

In our previous analysis (1; 2) we used an explicit parameter scan, which was limiting due to 
the large number of parameters. The new analysis uses the MultiNest software 2 (6), an advanced 
Bayesian package for multi-parameter fitting; among its advantages are that it does not require 
specification of a step-size, has handling of multi-modal posteriors, and allow computation of the 
statistics of functions of the parameters. We found it fast and reliable for the problem at hand. The 
ranges of the parameters have to be specified. As described below, we fit to eq.(|6. 1|), and there are 
a total of 10 parameters to be fitted. The parameters in this application have rather limited ranges 
of reasonable values, and a flat prior over a prescribed range is sufficient for the present purpose; 
more complex priors are straightforward to include if required. MultiNest outputs full parameter 
chains, mean and standard deviations of each parameter, and mean and standard deviations of user- 
defined functions of the parameters, here chosen to be the formula representing the CR spectrum 
as specified above. 

4. y-ray production functions 

Recent reviews of hadronic y-ray production can be found in (7; 1; 2; 4; 8). Here we use the 
QGSJET-II-04 model (9; 7). In addition to photon production in proton-proton collisions, we have 
to account for the contribution to the photon yield by proton-helium, helium-proton, and helium- 
helium interactions. Since QGSJET-II includes nuclei-nuclei interactions, the parametrization of 
(7) can be used directly for all four reaction channels above the transition energy. For A > 4 we 

! An alternative method, in principle preferable, would be to perform a parameter-free deconvolution of the CR 
spectra, with some regularizing prior; the present method is an intermediate approach. 

^available at http://ccpforge.cse.rl.ac.uk/gf/project/multinest 
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apply a nuclear enhancement factor following (8). For proton energies below 20 GeV we use the 
y-ray production functions described in (1; 2) ; these are for proton-proton collisions only, so we 
again apply the factor given in (8) for p-He, He-p, Fie-Fie and heavier nuclei. The helium abundance 
of interstellar gas is 0.1 by number, and the interstellar gas composition for heavier nuclei uses (10). 
CR helium is taken into account explicitly, and the CR composition for heavier nuclei is based on 
ACE and CREAM. The contribution from CR and ISM with A>4 is about 10% relative to A<4, 
and is dominated by CR, much less coming from heavy gas nuclei. 

Bremsstrahlung from CR leptons is computed using the corresponding GALPROP routine. 

5. Synchrotron and direct constraints on electrons 

Synchrotron emission from electrons and positrons in the interstellar magnetic field provides 
essential constraints on interstellar leptons, independent of solar modulation. For a survey of ex¬ 
perimental data and theoretical arguments see (11; 12; 13). We are not concerned here with the 
injection spectrum of electrons, or how the interstellar spectrum is affected by propagation, but 
just the observational information on the interstellar spectrum. As found in (11) the synchrotron 
brightness temperature spectral index /J changes from about 2.5 to 3 in the range from 100 MHz to 
several GHz. The relation of /3 to the electron index a is j3 = 2 ± (a — l)/2, so this corresponds to 
a steepening of the interstellar electron spectrum from index 2 to 3 at a few GeV electron energy. At 
high energies and frequencies, this is consistent with the Fermi-LAT direct measurements of elec¬ 
trons + positrons, which give an index 3.08 ±0.05 (14) and for electrons an index 3.19 ±0.07 (15), 
the electrons fully dominating at low energies where bremsstrahlung is important; this is also con¬ 
sistent with microwave synchrotron emission measured with Planck (16). Here solar modulation is 
negligible so that this consistency is a requirement, assuming the local measurement is typical of 
the interstellar medium near the solar position in the Galaxy. At low energies, only synchrotron is 
a reliable tracer of the interstellar electron spectrum, due to the very large solar modulation (factor 
or 10 or more). The relation of synchrotron to the electron spectrum is a complicated function (11), 
but it is useful to write a simplified version to illustrate the nature of the constraints on the energy 
of the electron spectral index break. We use the formula from sec 2.1.1 and 2.3 of (11) and the 
synchrotron data from that paper. We have for the synchrotron peak frequency for a electron en¬ 
ergy E: v max = 0.29v c = 0.29^^^ [^] 2 x 240 MHz. so E = y/ hm GeV ■ Note that 

E depends on \/v max and \[~B so is robust against uncertainties in these quantities. The observed 
break v is in the range 500-5000 MHz, B is between 5 and 10 pG, using a broad conservative 
range. Hence taking extremes of v/B, 2 < Ebreak < 10 GeV. A narrower range is also reasonable: 
v = 500 — 2000 MHz, B = 5 — 8 pG, giving 2.6 < Eb rea k < 6.6 GeV. The upper limit is consis¬ 
tent with Fermi-LAT electrons which show no break down to 7 GeV. However other experiments 
(AMS01, PAMELA) show the break must be below 3 GeV, since there is no break observed above 
this energy and solar modulation has the effect of enhancing the break in the directly measured 
spectrum. Hence the adopted break energy constraint from synchrotron and direct measurements is 
1-2.5 GeV. The synchrotron constraint is essential for the low end and the index below the break. 
Since we use a smooth function for the break, the relation between the break parameter and the 
actual spectrum is not exact. We use the Fermi -LAT electron spectrum as measured for the spec¬ 
trum above the break, including the Fermi-LAT normalization at 100 GeV and index range 3.1-3.2. 


4 






Local interstellar cosmic-ray spectra 


A. W. Strong 


Below the break, the experimental range of /3 is 2.4 to 2.6, so from the above formula we use the 
range a = 1.8 — 2.2 as the constraint from synchrotron. 

6. Parameterization of spectra 

As explained in (1; 2) we use a spectrum in the form of density per unit momentum n(p), 
since this is expected to have a power-law shape in the theory of diffusive shock acceleration, and 
hence serves as a suitable basis form on which propagation and other effects are superimposed. 
The local emissivities require a turnover in the proton spectrum towards low energies. A possible 
parameterization of CR spectra has a single sharp break in the spectral index, and this is often used; 
it was used in for example in (1; 2). However this is unphysical, and a smoother function is more 
appropriate. Hence smooth breaks in protons, helium and electrons are modelled; the functional 
form used is n(p) °< 1 / + {jp) ai ^ S ] 5 where are the indices below and above the 

break respectively for oq < Ch, and the spectrum breaks around a momentum centred on /?/„.; the 
parameter 8 controls the sharpness of the break: smaller 8 produces a sharper break, typical values 
are 8 =0.5 - 1.5. It converges to the given power-laws at low and high p, with a smooth transition. 
This is the same form as used for CR in supernova remnants in (17), written in a symmetrical form 
in the indices for clarity. This function has one more parameter (8) than for a perfectly sharp break; 
it simply makes the break smoother in a controlled way. Given a reference value n re f at p re f the 
normalized spectrum is: 

n{p) = n ref x [(^)“>/ 5 + (^)«*/«]« / [(_P_)«i I s + (6.1) 

Pbr Pbr Pbr Pbr 

The normalization n re f(p re f) is treated as a free parameter; there are thus 5 parameters each 
for protons and leptons, making 10 parameters in total. Since gamma rays cannot distinguish an 
origin in CR protons from CR helium, some assumptions have to be made to break the degeneracy. 
The ratio of protons and helium is fixed to that measured by PAMELA at 100 GeV, where solar 
modulation is absent while the experimental error is small. We take a reference proton momentum 
of 100 GeV, with a wide prior range since the normalization is to be determined from the gamma- 
ray data only. The lepton (electron + positron) spectrum is only marginally constrained by gamma 
rays, so we rely on direct measurements for the normalization and high-energy spectral index, 
allowing a range with uniform prior. We use the Fermi-LKT lepton flux at 20 GeV (14). For the 
low-energy (below about 1 GeV) spectral index we rely on synchrotron radiation constraints as 
described later. The prior range for the break position is based on both synchrotron and direct 
measurements. 

7. Results and Discussion 

We illustrate the results using the raw emissivities accounting for energy dispersion; using 
the corrected values without applying dispersion leads to very similar results. The fit results are 
given in Table 1 and Fig 1. The parameters are highly correlated so the individual values are just 
indicative, while the resulting spectrum error band (Fig 1) uses the full posterior in all parameters. 

The proton spectrum spectral index above the break is recovered in good agreement with direct 
measurements (beyond the effects of solar modulation); note that it has been determined from the 
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Table 1: Summary of model fits to equation 6.1. Entries are prior range, posterior mean and standard 
deviation. The proton parameters are constrained by the y-ray emissivities, while the lepton parameters re¬ 
flect mainly the prior from synchrotron and direct measurements. The parameters are highly correlated and 
degenerate, so the resulting spectrum derived from the full posterior (Fig 1) is preferred to the individual pa¬ 
rameters. The CR density n re f is multiplied by (c/4k) to give a flux in the usual units quoted in experiments. 
p re f= 10 5 MeV for protons, 2 x 10 4 MeV for leptons. 


Parameter 

range: min 

max 

mean 

std 

units 

Protons 

(c/4n)n ref 

1 x 1(T 9 

20 x10~ 9 

6.4 x 10~ 9 

0.3 x 10~ 9 

cm 2 sr 1 s 1 MeV 1 

a\ 

2.2 

2.7 

2.37 

0.09 


«2 

2.6 

3.5 

2.82 

0.05 


8 

0.05 

1.0 

0.5 

0.1 


Pbr 

1000 

10000 

5870 

2200 

MeV 

Leptons 

(c/4n)n ref 

1 x 10~ 9 

3 x 10~ 9 

2.2 x 10~ 9 

0.5 x 10~ 9 

cm 2 sr 1 s 1 MeV 1 

a\ 

1.8 

2.2 

2.0 

0.1 


«2 

3.1 

3.2 

3.15 

0.03 


8 

0.05 

1 

0.47 

0.25 


Pbr 

500 

2000 

1130 

4067 

MeV 


gamma-ray data alone. The spectral break is significant, the index being smaller by ~ 0.5 below the 
break; this is where solar modulation affects the direct measurements, as can clearly be seen in the 
difference between PAMELA, AMS01 and the interstellar spectrum. Note that we address a break 
in the momentum spectrum, not kinetic energy; the latter representation is not appropriate, see (1; 
2). The overall normalization of the proton spectrum is about 30% larger than direct measurements, 
even in the high-energy region free from modulation. This can be a real physical effect since there 
is no guarantee that the direct measurements represent the average in the local region within one 
kpc probed by the gamma rays. There are spectral variations between the local region and the 
rest of the Galaxy (18; 19), and the region probed by the enrissivity analysis might be not so 
“local” in this sense. However there are sufficient uncertainties in the analysis to make such an 
interpretation premature at this stage: hadronic cross-sections have still significant uncertainties 
especially for CR and target nuclei with A > 1. Uncertainties in the gas column densities including 
gas not traced by HI or CO are still significant, even though the emissivity derivation is for atomic 
hydrogen well traced by the 21-cnr line, correlations between the gas phases does not allow full 
separation of the components. Ionized hydrogen is accounted for in the analysis, but again cannot 
necessarily be fully separated from the other components. We note that (5) implicitly finds a similar 
trend in the normalization: that analysis included proton and helium direct measurements in the 
fitting so is more constrained to agree with those. In fact their proton spectrum is slightly above 
direct measurements, and their predicted emissivities are slightly below the measurements; so to 
reproduce the emissivities the proton spectrum would have to be increased further relative to direct 
measurements, in accord to what we find here. Different cross-sections and analysis techniques (see 
below) can account for any remaining differences. One important use of the interstellar spectrum is 
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Figure 1: Cosmic-ray and emissivity spectra derived from model fitting. Yellow band shows model range. 
Model ranges are 1 standard deviation on the parameterized synthetic spectra. Upper: Measured and derived 
cosmic-ray proton spectra. Data are AMS01 (asterisks) and PAMELA (diamonds). Spectra are multiplied 
by p 2 , and also times p 2 8 to better show the break. Lower left: Measured and derived cosmic-ray electron 
spectra. Data are AMS01 (asterisks), PAMELA (diamonds), and Fermi -LAT (squares). The CR density n(p) 
is multiplied by (c/An) to give a flux at high momenta in the usual units quoted in experiments. Lower right: 
Fermi-LAT emissivity data (vertical bars) and model, with red and green curves showing the hadronic and 
leptonic bremsstrahlung contributions; the yellow band shows the total. 


for solar modulation studies; in this case it would be appropriate to use the spectral shape, which is 
well determined, but renormalize the spectrum to agree with direct measurements at high energies. 
This would attribute all the normalization difference to the factors mentioned above. 

The interpretation of the curvature in the proton spectrum is beyond the present work, but it 
is generally consistent with expectation based on CR propagation implied by the B/C ratio, which 
peaks at a few GeV, implying either a break in the diffusion coefficient, effects of diffusive reaccel¬ 
eration, or convective transport. Any of these effects will cause a low-energy break in the primary 
spectra including that of protons, superimposed on a power-law injection spectrum (20). Note the 
apparent difference in index between the local proton spectrum and the large-scale Galaxy (18; 19); 
the latter spectrum is evidently harder. The origin of the difference is not clear at present, but it 
is interesting that the index from local gamma rays and local direct measurements are consistent. 
There are further suggestions of local CR spectral variations (21). 
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8. Comparison of analyses 

The present analysis differs from that of (5) (Casl5) in various ways; both approaches are valid 
and are complementary, emphasizing different aspects of the subject. The main differences can be 
summarized as follows: Casl5 performs an iterative correction for energy dispersion, we include it 
in the response to the raw emissivities. Casl5 includes direct measurements of protons and helium 
in the fit, combined with the emissivities. We use only y-ray data for the hadronic contribution, and 
adopt a He/p ratio from direct measurements, assuming equal spectral shapes for p and He. Casl5 
uses a velocity term to obtain the low-energy turnover in the proton spectrum; we use a smoothly 
broken power law for more flexibility. Casl5 uses a flux spectrum in kinetic energy, we use a 
density spectrum in momentum. Casl5 uses cross-sections from Kamae; we use a combination 
from Dermer and QGSJET-II. Casl5 uses a lepton spectrum from direct measurements combined 
with the y-ray emissivities, we combine direct measurements with constraints from synchrotron. 
Casl5 uses Minuit for the fitting, we use a Bayesian scheme with MultiNest. 
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